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Stable law distributions occur in the description of the linear dielectric behavior of polymers, the motion 
of carriers in semi-conductors, the statistical behavior of neurons, and many other phenomena. No accurate 
tables of these distributions or algorithms for estimating the parameters in these relaxation models exist. In this 
paper we present tables of the functions 
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together with related functional properties of zQ a (z). These are useful in the estimation of the parameters in 
relaxation models for polymers and related materials- Values of the integral QJz) are given for a=0.01, 
0.02(0.02)0.1(0.1)1.0(0.2)2.0 and those of V a (z) are given for a=Q.0(0.01)0.1(0. 1)2.0. A variety of methods was 
used to obtain six place accuracy. The tables can be used to sequentially estimate the three parameters appearing 
in the Williams-Watts model of relaxation. An illustration of this method applied to data in the literature is given. 



1. Introduction 

Stable law distributions and functionals of these distri- 
butions play an important role in a variety of scientific 
areas. They originated in a study of observation errors 
by Cauchy [l] 1 , and many of their mathematical proper- 
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ties were elucidated by Levy [2], and Khintchine and 
Levy [3j. Applications of stable laws have appeared in 
the context of models for the broadening of spectral 
lines [4], fluctuations in stock market prices [5], statisti- 
cal properties of neuronal activity [6], the motion of 
carriers through amorphous semiconductors [7], as well 
as in a variety of chemical physics problems [8]. Inte- 
grals of stable law densities have appeared in still an- 
other area of application: the theory of mechanical and 
electrical relaxation processes. The suggestion that re- 
laxation processes in some electrical and mechanical 
systems could be described in terms of stable distribu- 
tions can be traced back to experiments by Weber [9,10] 
and the Kohlrausches, father [11] and son [12]. Consid- 
erable recent interest in relaxation processes describable 
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in terms of stable laws has been inspired by experiments 
performed by Williams and Watts on polymers [13,14], 
Many other examples of such applications to polymers 
and glasses occur in the literature, cf., for example, the 
work of Moynihan, Boesch, and Laberge [15] and 
Bendler [16]. The marriage between the mathematics of 
stable processes and the theory of dielectric relaxation 
processes has recently been proposed by Schlesinger 
and Montroll [17] who used the methodology of con- 
tinuous time random walks [18] to try to model the 
underlying physical processes. 

The relaxation of linear systems can be characterized 
by a frequency dependent dielectric constant or the ap- 
propriate mechanical analogue. This constant is gener- 
ally expressible as 



€(ft))-6(oo) 

£(0)-e(oo) 



= e' a (o>)-i£'n («)=-J V*'^ A 



(1) 



where <£>(/) is the relaxation function which character- 
izes the particular material properties and e'„((oJ and 
e" (w) are the components of the normalized dielectric 
constant. It is found that for many polymers and glasses 
<J>(0 can, to a good approximation, be chosen to have the 
form 



4>(t)=exp[-(t/rT] 



(2) 



where r and a are constants that depend on the material. 
Data on polymers and glasses generally lead to a values 
in the range 0.2 to 0.8. Equations (1) and (2) together 
imply that, for z = atT, 

e' n (o))= 1 — z e~" a sm(zu)du = 1 — ttzVJz) 

e"(o>)=2 e~ ua cos(zu)du=TrzQJz) (3) 

Jo 

where QJz) and VJz) are the standard integrals 
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stable density they were able to check their results 
against Holt and Crow's tables [20] for a— 1/4, 1/2, and 
3/4. However, these tables are given to four places in 
fixed arithmetic. The variety of methods used in gener- 
ating the tables has resulted in a lack of uniform accu- 
racy. 

Because of the wide range of applications of the func- 
tions QJz) and VJz), and because there appear to be no 
accurate tabulations against which approximations can 
be checked, we present here tables of these functions for 
a=0.01,0.02(0.02)0.1(0.1)1.0(0.2)2.0 for QJz), and for 
a=0.0(0.01)0. 1(0. 1)2.0 for V a (z), accurate to six places 
in floating point. 



2. Numerical Methods 

Three methods were used to evaluate QJz) and VJz) 
to 10 significant figures, which were truncated and 
rounded to six figures for the present tables. The first is 
the evaluation of a convergent series representation, the 
second is the evaluation of a divergent series, and the 
third is numerical integration in the region in which the 
divergent series does not yield the required accuracy at 
all, and the convergent series requires use of an unrea- 
sonably large number of terms. For Q a (z) the following 
series converges for 0<ct<l 



e^)=^|(-ir'^^ S m(^f). (5 ) 



This series diverges when l<a<2, but it is an asymp- 
totic series, [21,22], and can be profitably used in this 
range for computation at sufficiently large values of z. 
The corresponding series that converges in l<a<2, 
but diverges when 0<a<l, is 



^=^2,<- 1 >• +, lb32j/ ^2( "-' , 



(2«-2)! 



(6) 



Similar series can be derived for V a {z). These are 



VJz)=- fV^shMzwyw. (4) 

The function Q a (z) can be identified as a representation 
of a stable law density. In a recent paper Montroll and 
Bendler have presented a number of approximations to 
the function QJz) for values of a useful for polymer 
applications [19]. Because QJz) can be identified with a 
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which converges for 0<a<l and diverges when 
l<a<2, and 



r (^) 



z 2 "- 1 (8) 
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which diverges in 0<a<l and coverges in l<a<2. 
When the appropriate series from eqs (5-8) requires a 
large number of terms to provide accurate results it is 
convenient to make the substitution u = tan0 in eq (4) to 
find the following integral representations: 



77 Jo 



71/2 rift 

e- tanafl cos(ztan0) aV 



cos 2 



(9) 



77 Jo 



~ Uaae sin(ztan0) ~^- 



cos 2 



which may be evaluated numerically. 

Checks on the numerical calculations can be made for 
a few values of a for which both QJz) and V a (z) can be 
expressed in terms of known functions. For a= 1/2 one 
finds 



Q m {z)=2Trp> {[i-C(p)] cos^) 
|-5(p)]sin(^)} 



+ 



(10) 



-1/2 



where 

P=(2ttz)- 

C(p)=j P cos(^fyt, S(p)=j\m^pjdt (11) 

in terms of Fresnel integrals. The conjugate function, 
V\n{z), can also be written in terms of p as 

V U2 (z)=2p 2 l-7rp{[i-5(p)]cos^) 

- J-C(p)]sin(^)}]. (12) 

The functions in curly brackets in the expressions for 
Q\a(z) and V\ n {z) are tabulated in Abramowitz and Ste- 
gun [23]. One can also verify the following special cases: 

gl(Z) = 7T(l+Z 2 ) ' W =2 &( Z > 



^(z)=(47r)- 1/2 exp(-z 2 /4) s 



Cz/l 

V 2 (z) = it- ' exp( -z 2 /4) extft 2 )dt 
Jo 

the expression for F 2 (z) being proportional A ~ 



to 



(13) 
Dawson's 



integral of argument z/2 [23]. In addition, Zolotarev 
[24] has found the value of Qmiz) in terms of Whittaker 
functions, but we have not evaluated his expression. 



3. Computational Notes 

All of the computations were performed in double 
precision in FORTRAN 77 on a Perkin-Elmer 3230 
computer. Each complete set of tables, Q a (z) and V a (z), 
took between 1 and 1.5 minutes of CPU time. Tables 1 
and 2 indicate when the appropriate series was used for 
computing Q a (z) and V a (z), respectively, as well as 
when numerical integration was required. Two points 
should be noted with respect to the two tables: 1) table 
entries for the critical values of z will change somewhat 
if other than six-digit accuracy is required, and 2) entries 
for critical z in tables 1 and 2 assume only restricted 
values, namely those values of z that appear in tables 3 
and 4 (positioned at the end of this paper because of their 
length), where the Q a (z) and V a {z) are tabulated. Note 
that for some values of a the critical z values produce a 
range of overlap (for example, for a = 1.9 in table 2); in 
these cases there is no need for numerical integration. 
But in other instances there is a gap between the critical 
values of z for use of the respective series, which may 
simply be due to the fact that no values of z are tabulated 
in the intermediate range. Hence, in these instances the 
calculation of QJz) or VJz) for intermediate values of 
Z may or may not require numerical integration. Since 
double precision was used in the calculations, some care 



Table 1. Critical values of z for the use of eqs (5) and (6) or 
numerical integration to evaluate Q a (z). 









Use 




Use 


Use 


numerical 


a 


eq(6) 


eq(5) 


integration 




for z < 


for z> 


for z = 


<0.3 


— 


0.001 


— 


0.4 


0.002 


0.004 


0.003 


0.5 


0.01 


0.02 


— 


0.6 


0.05 


0.10 


— 


0.7 


0.1 


0.15 


— 


0.8 


0.25 


0.35 


0.30 


0.9 


0.50 


0.65 


0.55,0.60 


1.1 


1.0 


2.0 


1.5 


1.2 


2.0 


2.5 


— 


1.3 


3.0 


3.5 


— 


1.4 


4.0 


4.0 


— 


1.5 


4.5 


5.0 


— 


1.6 


5.5 


5.5 


— 


1.7 


6.0 


6.5 


— 


1.8 


6.5 


7.5 


7.0 


1.9 


6.0 


8.5 


6.5-8.0 
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Table 2. Critical values of z for the use of eqs (7) and (8)or 
numerical integration to evaluate V a {z). 









Use 




Use 


Use 


numerical 


a 


eq(8) 


eq(7) 


integration 




for 2< 


for z> 


forz = 


<0.3 


— 


0.001 


— 


0.4 


0.001 


0.01 


0.002-0.005 


0.5 


0.01 


0.025 


0.015,0.02 


0.6 


0.045 


0.06 


0.05,0.055 


0.7 


0.1 


0.15 


— 


0.8 


0.25 


0.35 


0.3 


0.9 


0.5 


0.65 


0.55,0.6 


1.1 


1.0 


2.0 


1.5 


1.2 


2.0 


2.5 


— 


1.3 


2.5 


3.0 


— 


1.4 


3.5 


4.0 


— 


1.5 


5.0 


5.5 


— 


1.6 


5.5 


6.5 


— 


1.7 


6.5 


6.0 


— 


1.8 


7.5 


7.5 


— 


1.9 


8.0 


7.5 


— 


2.0 


8.5 


8.0 


— 



Watts relaxation function. Although it appears that all 
three parameters must be fit simultaneously, we will 
show that a small extension of tables 3 and 4 allows one 
to estimate the parameters separately. Since the most 
frequently measured property from which it is possible 
to determine the parameters a, r, and A = e(0) — 1 is the 
dielectric loss function, i.e., €"(&>), we restrict our dis- 
cussion to this function. As can be seen from eq (3), the 
dielectric loss function is proportional to g a (z)=zQ a (z). 
This function is unimodal as a function of z as illustrated 
in figure 1. Several parameters can be defined that char- 
acterize the shape and properties of gjz). These include 
z m {a), the value of z at which g a (z) is maximized, 



M(a) = max gjz) =gjz m (a)) 



(16) 



the height of the peak, and for 0< 1 two sets of abscissas, 
{zi(6,a)}, {z t (0,a)}. These are, respectively, the values 
of z on the leading and trailing edges of g a (z) at which 



was required in the calculation of gamma function. We 
used the formula 



1 2a 

fOH^) = 2 a ^ 



(14) 



n=0 



for x=y + 1<3.21, and 



r(;c)=(27r) 1/2 ;c 



1/2„x-1/2--jc 



X c "*' 



(15) 



for x > 3.21. The coefficients {a n } and {c„} appearing in 
these formulae are to be found in the book by Luke [25]. 
The accuracy of formula [14] increases as x— >0, while 
the accuracy of formula (15) increases as x— >oo. For 
x = 3.21, both formulas yield the same accuracy. 

When the series of eqs (5-8) were used to evaluate 
QJz) or V a (z), the number of terms used was a max- 
imum of 150, which occurred near the switchover re- 
gion for the two series. Most entries in the tables never 
required the evaluation of more than 20 terms. The nu- 
merical integrations of eq (9) were performed using 
Simpson's extended rule with a step size of 7r/40,000. 



4. Application to Polymer Physics 

The expressions in eq (3) allow us to estimate the three 
parameters that characterize polymers whose dielectric 
properties are described by eq (2) which is the Williams- 



gJz)=BM(a). (17) 

These parameters are all illustrated in figure 1. For later 




Figure 1-A typical curve of gjz) =zQ a (Z) together with 
distinguished points useful for parameter estimation. The point Z m 
is the maximum, and Z/(0,a) and Z t {Q,a) are the two solutions to 
g a (Z)=eM(a). In the figure 6 was taken equal to 1/2. 
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reference we give values of z m (a) and M(a) in table 5 for 
a =0.05(0.05) 1.00. Table 6 contains values of 



f(0,a)=\n[ Zl (e,a)/z m (a)\ 



(18) 



for a =0.05(0.05) 1.00 and 0=0.1(0.1)0.9 which will be 
used for data analysis, as will be explained in further 



Table 5. Values of z m (a) and the maximum M(a) 
for a=O.O5(0.O5)1.0O. 



z m (a) 



M(a) 



0.05 


0.5741 


0.9188(-2) 


0.10 


0.5885 


0.1833(-1) 


0.15 


0.6041 


0.2738(- 1) 


0.20 


0.6207 


0.3632(- 1) 


0.25 


0.6380 


0.4513(-1) 


0.30 


0.6561 


0.5380(-l) 


0.35 


0.6751 


0.6233(-l) 


0.40 


0.694S 


0.7071(-1) 


0.45 


0.7154 


0.7896(-l) 


0.50 


0.7370 


0.8705(-l) 


0.55 


0.7595 


0.9501(-1) 


0.60 


0.7830 


0.1028(0) 


0.65 


0.8075 


0.1105(0) 


0.70 


0.8329 


0.1180(0) 


0.75 


0.8592 


0.1253(0) 


0.80 


0.8863 


0.1325(0) 


0.85 


0.9142 


0.1394(0) 


0.90 


0.9425 


0.1462(0) 


0.95 


0.9712 


0.1528(0) 


1.00 


1.0000 


0.1592(0) 



detail below. The values of a cover the physically inter- 
esting range (0.2, 0.8) which has so far been found in 
polymers. 

The procedure for parameter estimation involves 
three successive steps. Since t, a, and A are initially 
unknown one can get an estimate of a independent of r 
and^ by first estimating the peak height, z m (d), from the 
data and then solving for a from eq (18). Since z=2tt/t, 
the parameter r drops out of eq (18) when the ratio of z's 
is taken. Although it appears that there are many inde- 
pendent estimates of a, each of which corresponds to a 
different 0, the data may not be equally useful at all 
values of 0. At high frequencies the underlying physical 
assumption that the system has a single degree of free- 
dom may be poor approximation. In particular, the data 
at values less than 0.5 may not fit the Williams- Watts 
model as well as data for #>0.5. When a satisfactory 
value of a is determined, using eq (18), we may find t 
from the relation 



T=^(d)/(27T/ ra ) 



(19) 



where f m is the frequency in Hz at which the peak max- 
imum occurs. Finally, since an estimate of the peak 
height will generally be available, the estimate of the 
parameter A can be found from 



A =max £ {2irf m r)/M{a) 



(20) 



where relevant values of M(a) are given in table 5. 



Table 6. Values of the function f(6,a) = ln[z,(6,a)/z m (a)]. 





a 


0=0.1 


6=0.2 


6=0.3 


6=0.4 


6=0.5 


6=0.6 


6=0.7 


6=0.% 


6=0.9 


0.05 


6.5303(1) 


5.0611(1) 


4.1593(1) 


3.4833(1) 


2.9237(1) 


2.4287(1) 


1.9654(1) 


1.5039(1) 


9.9471(0) 


0.10 


3.2670(1) 


2.5323(1) 


2.0813(1) 


1.7432(1) 


1.4633(1) 


1.2157(1) 


9.8394(0) 


7.5297(0) 


4.9815(0) 


0.15 


2.1793(1) 


1.6895(1) 


1.3888(1) 


1.1633(1) 


9.7667(0) 


8.1152(0) 


6.5691(0) 


5.0280(0) 


3.3273(0) 


0.20 


1.6353(1) 


1.2679(1) 


1.0424(1) 


8.7325(0) 


7.3322(0) 


6.0932(0) 


4.9332(0) 


3.7766(0) 


2.4999(0) 


0.25 


1.3085(1) 


1.0146(1) 


8.3420(0) 


6.9892(0) 


5.8690(0) 


4.8778(0) 


3.9497(0) 


3.0244(0) 


2.0025(0) 


0.30 


1.0902(1) 


8.4535(0) 


6.9504(0) 


5.8235(0) 


4.8905(0) 


4.0649(0) 


3.2919(0) 


2.5210(0) 


1.6697(0) 


0.35 


9.3382(0) 


7.2399(0) 


5.9523(0) 


4.9871(0) 


4.1881(0) 


3.4812(0) 


2.8195(0) 


2.1596(0) 


1.4308(0) 


0.40 


8.1600(0) 


6.3250(0) 


5.1993(0) 


4.3558(0) 


3.6578(0) 


3.0405(0) 


2.4627(0) 


1.8866(0) 


1.2502(0) 


0.45 


7.2385(0) 


5.6087(0) 


4.6095(0) 


3.8611(0) 


3.2421(0) 


2.6948(0) 


2.1827(0) 


1.6724(0) 


1.1086(0) 


0.50 


6.4962(0) 


5.0310(0) 


4.1333(0) 


3.4615(0) 


2.9061(0) 


2.4154(0) 


1.9565(0) 


1.4993(0) 


9.9427(- 1) 


0.55 


5.8838(0) 


4.5537(0) 


3.7397(0) 


3.1309(0) 


2.6282(0) 


2.1843(0) 


1.7694(0) 


1.3561(0) 


8.9980(- 1) 


0.60 


5.3686(0) 


4.1517(0) 


3.4078(0) 


2.8521(0) 


2.3937(0) 


1.9893(0) 


1.6116(0) 


1.2355(0) 


8.2024(- 1) 


0.65 


4.9280(0) 


3.8075(0) 


3.1235(0) 


2.6133(0) 


2.1919(0) 


1.8223(0) 


1.4765(0) 


1.1324(0) 


7.5236(- 1) 


0.70 


4.5459(0) 


3.5084(0) 


2.8764(0) 


2.4057(0) 


2.0183(0) 


1.6774(0) 


1.3594(0) 


1.0431(0) 


6.9377(- 1) 


0.75 


4.2106(0) 


3.2459(0) 


2.6594(0) 


2.2235(0) 


1.8653(0) 


1.5504(0) 


1.2570(0) 


9.6523(-l) 


6.4276(-l) 


0.80 


3.9134(0) 


3.0129(0) 


2.4669(0) 


2.0620(0) 


1.7298(0) 


1.4382(0) 


1.1667(0) 


8.9674(-l) 


5.9811(-1) 


0.85 


3.6477(0) 


2.8047(0) 


2.2951(0) 


1.9180(0) 


1.6093(0) 


1.3385(0) 


1.0867(0) 


8.3625(- 1) 


5.5887(-l) 


0.90 


3.4082(0) 


2.6170(0) 


2.1403(0) 


1.7886(0) 


1.5012(0) 


1.2496(0) 


1.0155(0) 


7.8269(-l) 


5.2431(-1) 


0.95 


3.1911(0) 


2.4470(0) 


2.0006(0) 


1.6721(0) 


1.4042(0) 


1.1699(0) 


9.5210(-1) 


7.3521(-1) 


4.9389(-l) 


1.00 


2.9932(0) 


2.2924(0) 


1.8738(0) 


1.5668(0) 


1.3170(0) 


1.0986(0) 


8.9559(-l) 


6.9315(-l) 


4.67l5(-l) 
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This estimation procedure was applied to dielectric 
loss data of Sasabe and Moynihan [26] on polyvinyl 
acetate at 66.7 °C. The relevant data are the following: 



1. 
2. 
3. 
4. 
5. 
6. 
7. 



logic/ 

2.701 
2.854 
3.033 
3.187 
3.560 
3.701 
3.857 



0.964 
0.980 
0.948 
0.884 
0.665 
0.587 
0.510 



0.901 
0.678 
0.599 
0.520 



The first three data points were used to estimate the 
peak location and height. These were found to be 



log, /« =2.839 



e"(max)=0.981 



following which the values of shown in data set were 
calculated. Equations (18-20) were then used to esti- 
mate a, r, and A for the last four data points. The re- 
sulting estimates are the following: 



4. 
5. 
6. 

7. 



a 

0.62 
0.60 
0.60 
0.60 



r(s) 



1.83(-4) 
1.81(-4) 
1.81(-4) 
1.81(-4) 



9.54 
9.54 
9.54 
9.54 



These results should be contrasted with the estimates by 
Sasabe and Moynihan, d=0.59 and r= 1.97X 10~ 4 s for 
the same set of data. The fitting procedure used by them 
required approximations suggested by Moynihan, 
Boesch, and Laberge [15], but these are of uncertain 
accuracy over the entire range of a because no accurate 
tables were available to check them. 
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